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An extremely useful evolution equation that allows systematically calculating the two-time corre- 
lation functions (CF's) of system operators for non-Markovian open (dissipative) quantum systems 
is derived. The derivation is based on perturbative quantum master equation approach, so non- 
Markovian open quantum system models that are not exactly solvable can use our derived evolution 
equation to easily obtain their two-time CF's of system operators, valid to second order in the 
system-environment interaction. Since the form and nature of the Hamiltonian are not specified in 
our derived evolution equation, our evolution equation is applicable for bosonic and/or fermionic 
environments and can be applied to a wide range of system-environment models with any factorized 
(separable) system-environment initial states (pure or mixed). When applied to a general model of 
a system coupled to a finite-temperature bosonic environment with a system coupling operator L in 
the system-environment interaction Hamiltonian, the resultant evolution equation is valid for both 
L = and L / 

cases, in contrast to those evolution equations valid only for L = L' case in 
the literature. The derived equation that generalizes the quantum regression theorem (QRT) to the 
non-Markovian case will have broad applications in many different branches of physics. We then give 
conditions on which the QRT holds in the weak system-environment coupling case, and apply the 
derived evolution equation to a problem of a two-level system (atom) coupled to a finite-temperature 
bosonic environment (electromagnetic fields) with L 7^ L^. 
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I. INTRODUCTION 

In many fields of modern sciences, one has to deal with 
open quantum systems in contact with their quantum 
surroundings or environments (reservoirs or baths) [H-Q . 
Most often, one concerns with only the system dynamics 
and the key quantity is the reduced system density ma- 
trix p{t) defined as the partial trace of the total system- 
plus-reservoir density operator pxit) over the reservoir 
degrees of freedom; i.e., p{t) — Trji[pT{t)]. The evolution 
equation of the reduced density matrix is governed by the 
reduced Liouville equation or called the quantum master 
equation that can be Markovian or non-Markovian. 

Two-time (multi-time) correlation functions (CF's) of 
an open quantum system are important physical quanti- 
ties. They can provide significant information about the 
system, whereas the single-time expectation values can 
not. For example, the two-time CF's of the electromag- 
netic field emitted by an atom are required for calculat- 
ing the fluorescence spectrum ilriM- The two-time CF's 
of the number of emitted photons give the information 
about the photon statistics and describe the behavior of 
photon bunching and anti- bunching 1-5]. The two-time 
CF's of the electric current through nanostructure de- 
vices are useful in the study of the transport properties 
of current fluctuations and noise spectrum [10, In 
the Markovian case, an extremely useful procedure to 
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calculate the two-time (multi-time) CF's for open (dis- 
sipative) quantum systems is the so-called quantum re- 
gression theorem (QRT) [H-Q that gives a direct relation 
between the time evolution equation of the single-time 
expectation values and that of their corresponding two- 
time (multi-time) CF's. So knowing the time evolution 
of the reduced density matrix of the system allows one 
to calculate all of the single-time expectation values and 
two-time (multi-time) CF's in the Markovian case. This 
is, however, no longer true in the non-Markovian case. 
For non-Markovian open (dissipative) quantum systems, 
the QRT is not vahd in general Although it 

is commendable to use the exact procedures [8l.[9l.ll4l [T5j 
to calculate directly the non-Markovian two-time CF's, 
not too many problems can be exactly worked out in this 
way. It is thus important that a procedure similar to 
the QRT can be developed and be applied to calculate 
the two-time CF's perturbatively for the non-Markovian 
open (dissipative) quantum systems. 

Recently, by using the stochastic Schrodinger equation 
approach and the Heisenberg equation of system oper- 
ator method, Alonso and de Vega derived (T8l - [20l | the 
evolution equations of the two-time (multi-time) CF's of 
the system operators for a general model of a quantum 
system coupled to a bosonic environment with a system 
coupling operator L in the system-environment interac- 
tion Hamiltonian. The evolution equations, valid to sec- 
ond order in system-environment coupling strength, were 
applied to calculate the emission spectra of a two-level 
atom placed in a structured non-Markovian environment, 
i.e, electromagnetic fields in a photonic band-gap mate- 
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rial [2l|. In the photonic band-gap material, the cor- 
relation function of the electromagnetic field (environ- 
ment) is highly non-Markovian, particularly within the 
edges of the bands, so the QRT is not valid even when 
the atom-environment interaction strength is weak. Al- 
though the evolution equations, derived in Refs. [Tsl - Eol l 
(Eq. (6) in Ref. [H, Eq. (31) in Ref. [19] and Eq. (60) 
in Ref. [13l), are very useful to calculate the time evolu- 
tion of the two-time (multi-time) CP's of the system ob- 
servables, the derivations of the evolution equations were 
presented only for an environment at zero temperature 
and only for a system state in an initial pure state. In 
Ref. jlS] , it was mentioned that it is possible to use the re- 
duced stochastic system propagator that corresponds to 
an initial state of the environment different from the vac- 
uum to evaluate the single-time expectation values and 
multi-time CP's with more general initial conditions. But 
they derived only a master equation that is conditioned 
on initial bath states and is capable of evaluating just the 
single-time expectation values of system observables for 
general initial conditions, both for an initial pure state 
and mixed state [l^. So strictly speakin g, t he two-time 
evolution equations derived in Refs. |18l420l | is applica- 
ble only for a zero-temperature environment. However, 
these equations were used to calculate the two-time CP's 
of system observables of dissipative spin-boson models 
at finite temperatures (l8l - l20j . This is possible only for 
the dissipative spin-boson models with Hermitian system 
coupling operators, L = (see also the discussions in 
subsection IIII Ap . 

In this paper, we use another commonly used tech- 
nique, the quantum master equation approach fT-^, to 
derive in the weak system-environment coupling limit an 
evolution equation of the two-time CF's of the system 
operators for non-Markovian open quantum systems in 
finite-temperature environments /or any factorized (sepa- 
rable) system- environment initial states (pure or mixed). 
This quantum master equation approach, different from 
those in Refs. [l8l - l20j . allows us to explicitly point out 
an important nonlocal environment (bath) memory term 
that vanishes in the Markovian case but makes the evo- 
lution equation deviate from the QRT in general cases. 
Since the form and nature of the Hamiltonian are not 
specified in our derived evolution equation, the evolution 
equation is applicable for bosonic and/or fermionic envi- 
ronments and can be applied to a wider range of system- 
environment models. When applied to a general model 
of a system coupled to a finite-temperature bosonic en- 
vironment, our resultant two-time evolution equation is 
valid for both Hermitian and non-Hermitian system cou- 
pling operator cases. This is in contrast to the two-time 
evolution equations derived in Refs. [l8l - [20j valid only 
for the finite-temperature bosonic environment case with 
a Hermitian system coupling operator L = L^. Our 
derived equation that generalizes the QRT to the non- 
Markovian case will have broad applications in many dif- 
ferent branches of physics. 

The paper is organized as follows. In Sec. [TTl we 



derive the non-Markovian finite-temperature evolution 
equations for one-time expectation values and two-time 
CF's for a general system-environment thermal model. 
We then apply the derived equations (valid to second or- 
der in the system-environment interaction) to a thermal 
bosonic environment and work out the explicit forms of 
the evolution equations in terms of system coupling op- 
erators in Sec. IIIIl The evolution equations for a non- 
Markovian bosonic environment has been presented in 
Ref. ^22] without any derivation. In this paper, the de- 
tailed derivation is given. We also discuss the connection 
of our derived two-time evolution equations with those 
presented in Refs. (T8l - [20j and give conditions on which 
the QRT may hold in the weak system-environment cou- 
pling case. In Sec. IIVI we apply the newly derived equa- 
tion to a problem of a two-level system (atom) cou- 
pled to a bosonic environment (electromagnetic fields), 
in which the system coupling operator is non-Hermitian, 
i.e., L ^ . A short conclusion is given in Sec. IVl 



II. DERIVATIONS OF EVOLUTION 
EQUATIONS 

A. Time-nonlocal and time-local quantum master 
equations 

Let us consider a general total Hamiltonian of the sys- 
tem and reservoir (environment) as 



H — Hs + Hi + Hji, 



(1) 



where Hs and H^ are the Hamiltonian for system and 
reservoir, respectively, and Hj is the interaction between 
the system and reservoir. It is convenient to go to the 
interaction picture in which 



dpT{t)/dt^-{i/h)[Hi{t),pT{t)i 



(2) 



where Hi (t) = e'"'>*/'^Hie- is exphcitly time- 
dependent, pT{t) = e'-^ot/hp^^^-^^-^iHot/h jg ^-j^g ^Q^g^j 

sity matrix operator at time t in the interaction picture 
and Hq — Hs + H^. Then the Liouville equation of 
the total density matrix pt(^) in the Schrodinger picture 
becomes 



dpT{t) 



dt ^ 
One can integrate Eq. ^ formally to obtain 



-iHot/h 



dprit) 



,iHat/h 



PT {t) = PT (0) 



dt' 







Hj{t'),pT (f) 



(3) 



(4) 



and substitute for px (t) inside the commutator in 
Eq. ([2]). The resultant equation is 



dpTjt) 
dt 



Hi{t),pT (0) 



Hi{t), Hi{t')rPT{t') . (5) 
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Equation ([5]) is still exact. Suppose initially 
Pt{0) = Pt{0) = p(0) (g) Rq, where Rq = 
exTp{-HB/kBT)/Tiji[exp{-HR/kBT)] is the initial ther- 
mal reservoir density operator. Then after tracing over 
the reservoir, Eq. ([5]) gives the master equation of the 
reduced density matrix in the interaction picture 



dm 

dt 



dt' 



Hi{t), Hi{t'),pT{t') ,(6) 



where we have considered a class of system-reservoir in- 
teraction Hamiltonian models for the thermal reservoir 
such that 



Trn[Hi{t)Ro\ = 



(7) 



to eliminate the first term in Eq. ([S]). 

The perturbative non-Markovian open quantum sys- 
tem theory may be categorized into two classes: the 
time-nonlocal (or time-convolution) 0, 23-29| and the 
time-local (or time-convolutionless) [a M. I25l43(i | meth- 
ods. This has been discussed extensively in the literature 
0, HE HI) HI] • The environment or reservoir by defini- 
tion is large and contains many degrees of freedom so 
that the influence of the system on the reservoir is small 
in the weak system-environment coupling case. As a con- 
sequense, to second order in system-environment interac- 
tion, the total density operator on the right hand side of 
Eq. ([6]) can be approximated to an uncorrelated (factor- 
ized) state as prit') = p{t') ^ Rq + 0{Hi) |3] since the 
products of two interaction Hamiltonians iJ/'s appear al- 
ready there. So in many textbooks [H-Il], the replacement 
of prit') with p{t') (g) i?o in Eq. ([6]) is performed under 
the so-called Born approximation. One then obtains 1,^ 



dm 

dt 



dt' 



Hi{t), Hi{t'),p{t')®Ro 



(8) 



Note that Eq. (|8]) is in a form of delayed integro- 
differential equation and is thus a time-nonlocal mas- 
ter equation. However, it can also be shown that an- 
other systematically perturbative non-Markovian master 
equation that is local in time [1, 0, [25l - [36| can be de- 
rived from the time-convolutionless projection operator 
formalism [3, HB-HM or from the iteration expansion 
method Q . Under the similar assumption of the factor- 
ized initial system-reservoir density matrix state and the 
use of Eq. (O , the second-order time-convolutionless mas- 
ter e quat ion in the interaction picture can be obtained as 

[i, 0,11113 



dm 

dt 



dt' 



Hi{t), Hi{t'),p{t)®Ro 



(9) 



We note here that obtaining the time-convolutionless 
non-Markovian master equation perturbatively up to 
only second order in the interaction Hamiltonian using 



the time-convolutionless projection operator technique 
0, HE H^ is equivalent to obtaining it by replacing p(f') 
with m in Eq. §1 ^U,^^^- One may be tempted to 
think that the second-order time-nonlocal master equa- 
tion ([8]) is more accurate than the second-order time-local 
(time-convolutionless) master equation ([9]) since besides 
the Born approximation, the (first) Markovian approx- 
imation of replacing p{t') with p{t) in Eq. ([5]) seems 
to be an additional approximation made on the time- 
local master equation. But it was shown that this may 
not be the case. In many examples 0, [25l - l27l . [29l . Isif . 
the time-convolutionless approach works better than the 
time-nonlocal approach when the exact dynamics is used 
to test the perturbative non-Markovian theory based 
on these two approaches. Here we will consider the 
second-order non-Markovian time-convolutionless (time- 
local) evolution equation in our derivation. 



B. Quantum regression procedure 

Quantum regression procedure or QRT indicates that 
for (Markovian) open quantum systems, the equations 
of motion (or evolution equations) for single-time expec- 
tation values of system operators are also the equations 
of motion for two-time (multi-time) CP's. Formally, the 
single-time or one-time expectation values of system op- 
erators can be written as 

{A (ti)) = Trs^fl, [A (ti) pT (0)] = Trsofl [A (0) pr (ti)] , 

(10) 

where A{ti) represents a general system Heisenberg oper- 
ator(s), priti) is the Schrodinger total density matrix op- 
erator at time ti, and the subscript oi S (E) R attached to 
the trace symbol of Tr indicates a trace over the Hilbert 
space of the total composite system The evolution 
equation of the one-time expectation values of system 
operators can then be obtained as 

d {A {ti)}/dti = Trs^R [A{dpT{ti)/dti)] (11) 
= Trs [A{dpiti)/dti)] , (12) 

where we have used the fact that A — A{0) is a pure 
system operator and Tr ii[dpT{ti)/dt] — dp{ti)/dt. For 
the two-time CF with ti > t2, one has 

{A (ii) B {t2)) = Tvs^R [A it,) B {t2) PT (0)] 

= Tts^r [Ae-'"*/^BpT {t2) e^^*/''](13) 

where t — ti—t2, and A — A{0) and B — B{0) are system 
operators. Let xt{0) = BpT{t2)- Then the two-time 
CF (USD becomes {A (ti) B (^2)) = Trs®/? [Axt {t)]. It is 
then equivalent to the expectation value of the operator 
A with respect to the effective density matrix operator 



XT{t) 



-iHt/h 



BpT (^2) e*^*/'' that satisfies the same 



Liouville evolution equation as prit) even though xt(0 
may not be a proper density matrix (i.e., positive-definite 
trace-conservative operator). The evolution equation of 
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the two-time CF can be formally written as 



obtain up to second order in the interaction Hamiltonian 



d{A ih)B{t2))/dh 



Trs<g,R[AidxTit)/dt)] (14) 
Trs [Aidx{t)/dt)] , (15) 



where the relations of the reduced operator x{t) — 
Trfl[xT(i)] and Ti R[dxT (t) /dt] = dx{t)/dt have been 
used. If the reduced master equations dx{t)/dt and 
dp{t) I dt had the same operator equation form, one might 
conclude that the structure and the form of the evolution 
equation of the two-time CF would be the same as those 
of the single-time evolution equation and thus the QRT 
would apply. In fact, it has been shown that the QRT 
is not valid in general [H, [l^ , but the QRT or regres- 
sion procedure is useful and correct for systems where the 
coupling to reservoirs is weak and the Markovian approxi- 
mation holds 0, [l^, [l3| . The main purpose of the present 
paper is to derive the non-Markovian finite-temperature 
evolution equation of the two-time system CF's using a 
quantum master equation approach, an approach differ- 
ent from those in Refs. [l^,[20|. Our equations, which are 
valid for both a Hermitian and a non-Hermitian system 
coupling operators and thus generalize the corresponding 
results in Refs. [I^II^Ij can be used to calculate the two- 
time CF's for any factorized (separable) system-reservoir 
initial state and for any arbitrary temperature as long 
as the approximation of the weak system-environment 
coupling still holds. 



C. Evolution equations in the weak 
system-environment coupling limit 

Let us proceed to first derive perturbatively the ex- 
plicit evolution equation of the single-time expectation 
values in the non-Markovian case. Here we consider the 
second-order non-Markovian time-convolutionless (time- 
local) evolution equation in our derivation. We wish to 
obtain an evolution equation, dpT{t) / dt, valid to second 
order in system-environment interaction Hamiltonian, to 
substitute into Eq. pT|) for the single-time expectation 
values and into Eq. for the two-time CF's. It is 
convenient to first go to the interaction picture and ob- 
tain a time-local (time-convolutionless) evolution equa- 
tion of the density matrix valid to that order. This can 
be achieved by the substitution of prit') — Prit) in the 
second term on the right hand side of the equal sign of 
Eq. © [i,0,[2i-[33|. To go back to the Schrodinger pic- 
ture, we substitute the resultant second-order equation 
obtained from Eq. ([5]) into Eq. ^ to obtain the evolu- 
tion equation, dpT{t)/dt. By substituting this equation 
dpT{t)/dt valid to second order in the interaction Hamil- 
tonian into Eq. (jlip . the evolution equation of the single- 
time expectation values then consists of three terms. The 
second term involves the first term on the right hand 
side of the equal sign of Eq. JS]), and will vanish on 
the conditions that pt(0) = Pt{0) — p{0) <8) Ro and 
TrR[Hi{t)Ro] = [Eq. ^] are satisfied. As a result, we 



d{A jti)) I 

= -^Trsgifl ([i7s, A]pT(ii)) 



1 



drTis^R 



{Hi{T-ti)[A,Hi]pT{ti) 

+ [Hj,A]Hj{T^ti)pT{ti) 
{t/h)TTs^R{{[Hs,A]}{t,)pTm 



S(g)R 



{Hj{T-t,)[A,Hj]}{t,)pT{0) 
^{[Hi,A]Hi{t -t,)}{ti)pT{0)] X^G) 



where we have transformed from the Schrodinger picture 
to the Heisenberg picture in the second equal sign and 
{AB}{t) = exp{iHt/h)ABexp{-iHt/h). 

Since XTit) and pxit) obey the same equations of 
Eqs. ([3]) and ((5]), at first sight, one may think that the 
two-time evolution equations, Eqs. (|14l) and (jl5p . are sim- 
ilar to the single-time evolution equations, Eqs. (jlll) and 
(fT2|l . and thus might be tempted to conclude that they 
have the same form of the evolution equations. Indeed, 
by using Eqs. p^ . ^ and (O, the first and third terms 
of the resultant equation derived from Eq. (fT4)) are similar 
to the right-hand side of the single-time evolution equa- 
tion (fTH]) with the replacement of pt{0) — > xt(— ^2) = 
-S(^2)pt(0) and with the change of the integration region 
from [0,ti] to [t2,ti]. Then we obtain 



-Trs^Ri{[Hs,A]}{t,)Bit2)pTm 
1 



drTrscsR 
{Hi{r - ti)[A,Hi]}{ti)B{t2)pTiO) 
^{[Hi,A]Hi{T~ti)}{ti)B{t2)pTiO) 



(17) 



However, a significant difference is that the expectation 
values for the second term does not vanish, i.e.. 



-i/h)TTs^R 



Ae 



-iHat/h 



Hiit),XTiO) 



(18) 



in the non-Markovian case, where t = ti — <2 in Eq. ([T; 
The reason can be understood as follows. The interaction 
Hamiltonian Hj(ti —t2) in Eq. (jlSp involves the environ- 
ment operators in the time interval from t2 to ti , and the 
effective density matrix operator xt(0) can be written 
as xt(0) = xt(0) ^BpTit2) = BUit2,0)pT{0)U^t2,0), 
where [7(^2,0) — e"^^^^/'^ is the Heisenberg evolution 
operator of the total Hamiltonian from time to t2. 
If the environment is Markovian where the environment 
operator CF's at two different times are ^-correlated in 
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time, then we may regard that the environment oper- 
ators in Hj{ti — are not correlated with those in 
U{t2,0). So the trace over the environment degrees of 
freedom for operator Hi{ti — and operator U{t2,Q) 
can be performed independently or separately. The trace 
of pT{t2) = U{t2,0)pT{0)U^t2,0) over the environment 
degrees of freedom yields the reduced density matrix 
p{t2) — Trfl[pT(^2)], but the trace of Hi{ti — ^2) van- 
ishes, i.e., Trfl[7?/(ti — t2)Ro] ~ 0, because of Eq. (O. 
Thus Eq. (fT5|) vanishes in the Markovian limit. But the 
situation differs for a non-Markovian environment as the 
environment operator in Hj{ti — ^2) niay, in general, be 
correlated with that in [/(t2,0). Therefore, the evolu- 
tion from pt(0) to pT{t2) under the influence of interac- 
tion Hamiltonian in the presence of the reservoir needs 
to be taken into account before the trace over the en- 
vironment is performed in Eq. ^TE\i. We emphasize here 
that it is this nonlocal environment (bath) memory term, 
Eq. (|18l) , that vanishes in the Markovian case but makes 
the evolution equation of the two-time CF's of the sys- 
tem operators deviate from the QRT. As we aim to ob- 
tain an evolution equation of the two-time CF's of the 
system operators, valid up to second order in the inter- 
action Hamiltonian, we need to find pT{t2) only up to 
first order in the interaction Hamiltonian. So substitut- 
ing pT{t2) = e"*-^"*2/ri^^(^2)giffot2/?i ^itji the expression 

PT {t2) ^Pt{0)-j^J dT [Hj (t) , pT {t2)\ (19) 

for xt(0) = Xt(0) = BpT{t2) in Eq. p^ . we then obtain 
up to second order in the interaction Hamiltonian (in the 
system-environment coupling strength) 



j-2l drT, 



S®R 



B e~'"°'^/^[Hi{T),pT{t2)V"°'^''' 



JHot/h 



= - y dr Tis<sB. {A [Hi , 

e-'""'/'' B[Hi{t -t2), PTit2)] e'""'^" 



(20) 



where the first order term in the interaction Hamiltonian 
coming from pt(0) term in Eq. (jl9p has been dropped 
because of Eq. ([7|) . Since the product of two Hj already 
appear explicitly in Eq. (j20p . we may then transform 
Eq. ([20]) into Heisenberg representation with the evolu- 
tion equation exp{iHt/h) w exp{iHot/h) + 0{Hj). Fur- 
thermore, by writing out the commutators explicitly and 
rearranging the Heisenberg operator terms, the resultant 
equation from Eq. (PH)) then becomes 



S®R 



(^{Hi{T^h)[Hi,A]}ih)B{t2)pT{0) 



-{[A,Hi]}{h){BHi{T - t2mt2)pT{0) 



(21) 



The first term in Eq. (j21[) is ready to combine with the 
second term in Eq. (|17p to extend the integration from 



to ti. Similarly, one may rewrite the last term in Eq. ((2T|) 
using the relation Bfl^/(r-i2) = Hi{T-t2)B+[B, Hj{t- 
12)] so that the first new term can be combined with last 
term in Eq. ((T7)) to extend the integration from to ii . 

Putting all the resultant terms together, we obtain the 
evolution equation of the two-time CF's valid to second 
order in the interaction Hamiltonian as 

d{A{h)Bit2))/dti 
= ii/h)Trs^R {{[Hs,A]}(ti)B{t2)pTm 



1 



drTr 



S^R 



fi' Jo 

(^{Hj{T~h)[A,Hj]}{h)B{t2)pTiO) 

+ {[Hi,A]HiiT - ti)}{h)B{t2)pT{0) 
1 rt2 



Jo 



drTr 



Si»R 



{ [Hi, A]}{h){[B, Hj{t t2)]}(t2)PT(0) .(22) 



Compared to Eq. ([16]), it is the existence of the last term 
in Eq. 1^ that makes the QRT invalid. Equation ((^ 
is the main result of this paper. The derivation is based 
on perturbative quantum master equation approach, so 
non-Markovian open quantum system models that are 
not exactly solvable can use our derived evolution equa- 
tion to obtain the time evolutions of their two-time CF's 
of system operators, valid to second order in the system- 
environment interaction. In the derivation of Eqs. (jl6p 
and (1^^ . we have also used the assumption of a factor- 
ized initial system-bath state pr(0) = Pt{0) = p(0) Ro 
and the condition of Ti fj[H j (t) Rq] = 0, Eq. ([7]), to elim- 
inate the first-order term. Since the form and nature 
of the Hamiltonians are not specified, Eq. can be 
used to calculate the two-time CF's for non-Markovian 
open quantum systems with multi-level discrete or con- 
tinuous Hilbert spaces, interacting with bosonic and/or 
fermionic environments. The procedure and the degrees 
of difficulty to apply Eq. (|22|) to a open quantum sys- 
tem model (by taking into account nonlocal bath mem- 
ory effects and tracing out the bath degrees of freedom 
for factorized system-bath initial states) to obtain the 
two-time CF's of system operator are similar to those for 
the evaluation of the reduced density matrix of a second- 
order time-convolutionless non-Markovian quantum mas- 
ter equation [e.g., Eq. ([9])]. We will explicitly apply the 
evolution equation ([22) to a general model of a quantum 
system coupled to a finite-temperature bosonic environ- 
ment in Sec. IHII and a specific model of two-level system 
in Sec. IIVI Open quantum systems coupled to fermionic 
reservoirs (environments) could, for example, be quan- 
tum dots or other nanostructure systems coupled (con- 
nected) to nonequilibrium electron reservoirs (electrodes 
or leads) in the electron transport problems [Sll, |38| - |45| . 
The evolution equation ([22) can also be used to calcu- 
late the non-Markovian two-time CF's in such systems. 
In summary, our evolution equation (j22l) can be applied 
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to a wide range of system-environment models with any 
factorized (separable) system-environment initial states 
(pure or mixed). 



III. EVOLUTION EQUATIONS FOR THERMAL 
BOSONIC BATH MODELS 

To make contact with Refs. [l8l - l20| . we consider a 
quantum system coupled to a bosonic environment with 
a general Hamiltonian of the form 

H = Hs + Y,^9x (L^ax + ^4) + H ?iWAaIaA(23) 

A A 

where the system coupling operator L acts on the Hilbert 
space of the system, a\ and a\ are the annihilation and 
creation operators on the bosonic environment Hilbert 
space, and g\ and uj\ are respectively the coupling 
strength and the frequency of the Ath environmental os- 
cillator. 

Applying Eq. to Eqs. (HH) and and after trac- 
ing over the environmental degrees of freedom for factor- 
ized (separable) system-bath initial states, we arrive at 
the second-order evolution equations of the single-time 
expectation values 

d{A{h))/dti 
= (zA)Trs {{[Hs,A]}{t,)pm 

+ / drTrs 
Jo 

(a*(ii - T){DiT - h)[A, L]}{ti)p{0) 

+a{h - t){[L\ A]L{t - h)}{h)p{0) 

+P*{h^T){L{r-h)[A,L^]}{h)p{0) 

+ f3{h - T){[i, A]D{r - h)}{h)p{0)) , (24) 

and of the two-time CF's 

d{A {h)B{h))/dti 
^{ilh)i:Ts{{[Hs.A]}{h)B{t2)pm 

+ [ drTrs 
Jo 

(^a*{ti~T){D{T~h)[A,L]}{h)B{t2)p{0) 

+a{h ~ T){[L\A]LiT - h)}{h)B{h)p{0) 
+I3*ih - t){L{t - h)[A, L^]}{h)Bih)piO) 

+ P{h - r){[L, A]L\t - t,)}{t,)B{t2)p{Q)) 

+ / drTrs 
Jo 

(a{t,^T){[L\A]}{t,){[BMr-m{h)p{Q) 
+ P{h - T){[L,A]}{t,){[B, D{t - hmh)p{())) (25) 



Here L{t) = exp{iHst/h)Lexp{—iHst/h) is the system 
operator in the interaction picture with respect to Hs, 
and 

a{h-r)=J2 ("A + 1) IffAl' e-*"^(*^-^), (26) 

A 

/3(*i-^)=E^A|5A|'e*'^^(*^-^) (27) 

A 

are the environment CF's with a{ti — r) = 
{'E\9\ax{ti)J2\' 9yai,iT)'^ ^ and /3{ti - t) = 

{j2\9\al{ti)J2x' 9ya\'{T)'j ^, where ixih) = 

oag"*'^^*^ and a\{ti) = a^e*'^^*^ are the environ- 
ment operators in the interaction picture, and the 
symbol {■■■)r denotes taking a trace with respect to 
the density matrix of the thermal bosonic reservoir 
(environment). The thermal mean occupation number 
n\ of the bosonic environment oscillators in Eqs. (j26p 
and (EH) is nA = (e'*"^/'=«^ - 

The evolution equations (|24l) and ([25)) for a non- 
Markovian bosonic environment have been presented 
in Ref. (22| without any derivation. In this paper, 
the detailed derivation of the evolution equations is 
given. Furthermore, the two-time CF evolution equation, 
Eq. p2|) . applicable for both bosonic and fermionic envi- 
ronments and applicable for more general form of system- 
environment interaction Hamiltonian has not been pub- 
lished in the literature yet. 

As mentioned, the two-time evolution equations de- 
rived in Refs. |18l - [20| is, strictly speaking, applicable 
only for a zero-temperature environment. However, these 
equations were used to calculate the two-time CF's of 
system observables of dissipative spin-boson models at 
finite temperatures. This is possible only for the dissipa- 
tive spin-boson models with Hermitian system coupling 
operators, L = . We will discuss this point in details 
in subsection IHI Al In contrast, our bosonic evolution 
equations, Eqs. ([M)) and (I25|) . are valid for both a Her- 
mitian and a non-Hermitian system coupling operators 
and can be used to calculate the two-time CF's for any 
factorized (separable) system-reservoir initial state and 
for any arbitrary temperature as long as the assumption 
of the weak system-environment coupling still holds. 

In Ref. [Hi, we used Eqs. (HD) and ^ to calculate 
the finite-temperature single-time expectation values and 
two-time CF's for a non-Markovian pure-dephasing spin- 
boson model of 

Hs = {nus/2)a,, L = n, = LK (28) 

Since the non-Markovian dynamics of this exactly solv- 
able pure-dephasing model can be cast into a time-local, 
convolutionless form and \L,Hs\ = 0, the results ob- 
tained by our second-order evolution equations turn out 
to be exactly the same as the exact results obtained 
by the direct operator evaluation. However, these re- 
sults significantly differ from the non-Markovian two- 
time CF's obtained by wrongly directly applying the 
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quantum regression theorem (QRT). This demonstrates 
the validity of the evolution equations (IM)) and (P5|) . But 
the system coupling operators L of this pure depha sing 
model [23] and the examples calculated in Refs. (18142 Clj | 
are all Hermitian, i.e., = L. So we will present in 
Sec. IIVI the calculations of one-time averages and two- 
time CF's for a thermal spin-boson model with L ^ L\ 
for which only the evolution equations ([M)) and (1^51) , 
rather than those in Refs. [18l - l2(j |. are applicable. 



d_ 

dt 



A. Comparison and discussion 

We discuss in the following the connection of our de- 
rived two-time evolution equation (|25|) with those pre- 
sented in Refs. fT8'-'2^. In Ref. [l^, a master equation 
conditioned on initial coherent states of the environment 
in Bargmann representation, {zq^z'q), was derived in the 
weak system-environment coupling limit. Provided that 
the whole set of the initial conditions of the system of 
interest, \iIj{zq)), and the statistical probability J{zo, Zq) where 
for the member |i/)(zg))|zo) of the statistical ensemble 
are known, this master equation with Zg = zq is capa- 
ble of evaluating the evolution of single-time expecta- 
tion values for general initial conditions, including ini- 
tially correlated mixed states between the system and 
environment. However, the evolution equations of two- 
time (multi-time) CF's of system observables, Eq. (6) in 
Ref. 18], Eq. (31) in Ref. jll and Eq. (60) in Ref. [20|, 
were derived for an initial vacuum state of the environ- 
ment and an initial pure state of the system of interest. 
As a result, these two-time evolution equations are valid 
only for a zero-temperature environment (if the system 
coupling operator is not Hermitian, i.e., L ^ see dis- 
cussions below) . Compared with the corresponding zero- 
temperature two-time (multi-time) evolution equations 
derived in Refs. [18-20], our finite-temperature two-time 
evolution equation (|25p is valid for any initial factorized 
(separable) states (pure or mixed) at finite temperatures 
and for both Hermitian and non- Hermitian system cou- 
pling operators. The extra terms containing the bath CF 
j3{ti — t) or P*{ti — t) are due to the finite-temperature 
thermal environment. If we take the zero-temperature 
limit by letting n\ — Q and thus f3{ti — t) — fi*{tx — t) — 
0, as well as consider the initial pure-state case by letting 
Trs[---p(0)] ^ (V(0)]---]V(0)), then Eqs.dMl) and ^ 
reduce exactly to their corresponding zero-temperature 
pure-state evolution equations in Refs. [l8l - l20j . 

However, calculations of the two-time CF's of sys- 
tem observables of dissipative spin-boson models in 
finite-temperature thermal baths (rather than zero- 
temperature baths) were presented in Refs. [l8l - l20| even 
though in their derivations of the two-time (multi-time) 
evolution equations, the bath CF is given in its zero- 
temperature form. 



This is only possible due to the reason that the system 
coupling operator is Hermitian, L = L\ in the ther- 
mal bath examples presented in Refs. [18l - l20l ]. One may 
understand this as follows. It is known that the finite- 
temperature density matrix operator of a thermal bath 
can be canonically mapped onto the zero-temperature 
density operator (the vacuum) of a larger (hypothetical) 
environment ^46, ,47] . By mapping the total Hamiltonian 
Eq. (|23p and an initial thermal state to an extended total 
system with a vacuum state, the finite-temperature prob- 
lem can be reduced to a zero-temperature problem, and 
the resultant pure state ■0* = \''l't{z* ,w*)) for the system 
of interest satisfies the following linear finite-temperature 
non-Markovian stochastic Schrodinger equation with two 
independent noises zl and vol [SIS]: 



—iHsijjt + Lz^Tpt ~ / dTa{t — r) 
Jo 

+L^w;^t-Lj\Tl3{t-T)^, 
z: = -zJ2Vn^glzle^-^\ 



/■= * * —iujxt 



(30) 



(31) 



(32) 



are two independent, colored, complex Gaussian noises 
with zero mean and satisfy 

M[ztZr\ = M[z;z*] = 0, M[z;zr\ = a{t - t);(33) 
M [wtWr] = M JWt <] =0, M [wlwr] = I3{t - t)(.34) 

Here, and w\ are coherent state variables of the ex- 
tended environment in Bargmann representation, A^J- • • ] 
denotes the statistical average over the Gaussian pro- 
cesses Zj* and Wf*, and the bath CF's a[t — T) and j3{t — T) 
are defined in Eqs. ([26|) and ([ST)) , respectively. In the 
zero-temperature limit T — > 0, the mean thermal occu- 
pation number of quanta in mode A approaches zero, i.e., 
n\ — )• 0. and thus (3{t — r) — )• and a{t — r) — >• ao(t — r) 
that is defined in Eq. (P^. In this case, the noises, 
Eqs. dsn) and become z^ = Ea fl^e*'^^* ^^d 

ml — Q, and the finite-temperature equation (I30p reduces 



(35) 



that has a zero mean and satisfies 




,2 - 

9x\ e 



JA(tl--r) 



(29) 



to the simple zero-temperature equation 

d Sip 
—ijjt = -iHs^t + Lz^ipt ^ y dra^it - t)j^. 

Now consider the case of a Hermitian system coupling 
operator L — The finite-temperature equation pop 
can be simplified considerably by introducing the sum 
process = z. 

M[utUr] =X1«] = 
A4[ufUr] — acs{t — t) 

^ a{t -t)+ Pit - t) 

cos[uj\{t - t)] 

X 

-isin[LOx{t-T)]}. (37) 



(36) 
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In terms of the single noise process u^, the hnear finite- 
temperature non-Markovian stochastic Schrodinger 
equation (|30|) for the case of a Hermitian couphng opera- 
tor L = takes the simple form of the zero-temperature 
equation (1351) with the replacement of by Uj and 
ao{t — r) by acs{t — r) = a{t — r) + f3(t — r) defined in 
Eq. (|37p . It is for this reason of the Hermitian coupling 
operator L = L'^ — ax m the dissipative spin-boson 
model with a thermal environment that the two-time 
CF's of the system observables can be evaluated with 
the evolution equations derived in Refs. [is 20]. In 
other words, if the system operator coupled to the 
environment is not Hermitian L ^ L\ the two-time 
(multi-tirne) differential evolution equations presented 
in Refs. |18l - l20l | are valid only for a zero-temperature 
environment. 

In contrast, our two-time evolution equation psp is 
valid for more general finite-temperature cases where 
the system coupling operator is not a Hermitian op- 
erator, i.e., L ^ L'^ . In the case of L 7^ L\ our 
two-time evolution equation contains additional finite- 
temperature j3{ti — t) and /3*(ti — t) terms which can 
not be combined and simplified to a simpler form as the 
zero-temperature evolution equation derived in Refs. (isl - 
[20j . For a Hermitian coupling operator L ~ L'^ , one can 
see that besides the replacement of a more general sys- 
tem state with an initial pure system state by letting 
(■0(0)1 ••• IV-CO)) Trs[- • •p(0)], the finite-temperature 
evolution equation (|25l) reduces to its zero-temperature 
counterpart in Refs. [18 20] with the effective bath CF 
given by acfi — (3{ti — t) defined in Eq. (|37|. 

This demonstrates explicitly why the zero-temperature 
two-time evolution equations derived in Refs. jl8l - l20| can 
be used to calculate the system operator CF's for a ther- 
mal spin-boson model with a Hermitian system coupling 
operator. 



B. Conditions for the QRT to hold 



As mentioned in subsection IIIBl QRT is a very use- 
ful procedure that enables one (in certain circumstances) 
to calculate for two-time (multi-time) CF's of system 
operators from the knowledge of the evolution equa- 
tions of single-time expectation values. One can notice 
that if the last two terms in Eq. ([25]) vanish, then the 
non-Markovian single-time and two-time evolution equa- 
tions (1^^ and (psj) will have the same form with the 
same evolution coefficients and thus the QRT can ap- 
ply. As expected, these two terms vanish in the Marko- 
vian case since the time integration of the corresponding 
(5-correlated reservoir CF's, a{ti — r) oc 5{ti — r) and 
(3{ti — t) (X S{ti — t), over the variable r in the domain 
from to t2 is zero as ii > ^2- So from Eq. (^5)) . in the 
weak system-environment coupling case the QRT holds 
when (i) [L\A] = or [B,L{t - ^2)] = at the zero 
temperature, (ii) at finite temperatures, in addition to 
condition (i), the following condition also needs to be 



satisfied: [L,A] = or [B,D{t - ^2)] = 0, (in) in the 
Markovian case where the bath CF's are (5-correlated in 
time. 

Note that in some models, certain CF's, which formally 
obey the QRT but with evolution equations coupled with 
those of other CF's that do not obey the QRT, may yield 
solutions different from those given by the QRT |19|]. 



IV. APPLICATION TO A THERMAL SPIN 
BOSON MODEL WITH i / 

To illustrate the usage of the equations we have de- 
rived, we apply them to the problem of a two-level sys- 
tem coupled to a thermal reservoir, in which L ^ . We 
consider Hamiltonian, Eq. (P^ . with Hs — {hhJA/'2)czi a 
coupling operator L — and a system-environment in- 
teraction Hamiltonian whose magnitude is small enough 
to be considered as a perturbation. Since the coupling 
operator L ^ \s not Hermitian, the two-time (multi- 
time) evolution equations derived in Refs. [T^-|20| are not 
applicable and our evolution equation (j25p should be em- 
ployed to obtain the time evolutions of the two-time CF's. 



A. Single-time expectation values 

Before calculating the two-time CF's, it is instruc- 
tive to obtain the master equation of the reduced sys- 
tem density matrix for the model. Transferring from the 
interaction picture back to the Schrodinger picture for 
Eq. © and using the general Hamiltonian, Eq. (|23l) . 
we obtain the second-order time-convolutionless non- 
Markovian master equation for the reduced density ma- 
trix p{t) as 

dT{a{t - t)[L'^L{t - t)p{t) ~ L{t - t)p{t)L'^\ 

+/3{t - t)[LL^t - t)p{t) -D{t- t)p{t)L] 
-fh.c.}, (38) 



where a{t — r) and /3(i — r) are defined in Eqs. (f26| and 
(|27p respectively, h.c. indicates the hermitian conjugate 
of previous terms, and an operator with a tilde on the 
top indicates that it is an operator in the interaction 
picture. The only real assumption used to obtain Eq. (j38p 
valid to second order in the system-environment coupling 
strength is the total density matrix factorizing at the 
initial time t — Q, pt{0) — p{0) <8) Rq- Taking L = a^, 
Lt cr+, then L{t) = CT-e"'"-**, and D{t) = a+e^'^^\ 
and substituting them into Eq. (j38p . we obtain 

dt 2 ^ '^^ 

-{Tiit) {a+a^p{t)-n_p{t)a+) 

+T2{t) {a^a+p{t) - a+p{t)a^) + h.c.} , (39) 
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where 

Ti{t)= f dTa{t - T)e+"^^(^'-^\ (40) 
T2{t)= [ dTl3{t-T)e'"^^'''-^\ (41) 

"'0 

The master equation p9p is a time-local and con- 
volutionless differential equation. The effect of the 
non-Markovian environment on the second-order mas- 
ter equation p9l) is taken into account by the time- 
dependent coefficients ri(ti) and T2{ti) defined in 
Eqs. (|40p and (|4ip instead of memory integrals. Like- 
wise, the evolution equations of single-time expectation 
values and two-time CF's of system operators valid to 
the second order are also expected to be convolutionless. 
Taking again L = o-_, = cr+, then L{t) = cr_e~*'^^*, 
and (t) = cr+e*"-**, and using the commutation relation 
between the Pauli matrices, we obtain straightforwardly 
from Eq. the following evolution equations of the 
single-time expectation values of system operators as 

d{a+{ti))/dti = iujA{(J+{ti)) 

-[rj(ti) + r2(ti)](a+(ti)), (42) 
d{a^{ti))/dti = -iujA{o--{ti)) 

-[T,ih)+T*{h)]{a^it,)), (43) 

d(a,(<i))Mi = - [Tiih) + r*(<i) + r;(<i) + r2(ti)] 

x(^.(ii)) 

-[T,ih)+Tlih)-T*ih)-r2ih)]. 

(44) 

Equations (|42p - (P^ can also be obtained directly from 
the master equation ((39)) through d{A{ti))/dti = 
Trs [A{dp{h)/dh)]. 

B. Two-time CF's 

For the evaluations of the two-time CF's of the system 
observables, we consider the following four cases. 

Case 1: [L\ A] = or [B,L{t)] = 0; and [L,A] ^ 
or [B,U{t)\ = 0. In this case, let A = cji^ B = 
Ui with i — +,— . Applying the commutation rela- 
tions between the Pauli matrices and the definition of 
Trs[(Tj(ti)(T^(t2)/o(0)] = {(Ji{ti)ai{t2)) to the right-hand 
side of the two-time evolution equation (^5]) , we then ob- 
tain 

d{(T+{ti)a+{t2)) / dti = iujA{(7+{ti)<T+[t2)) 

-[Tl{t,)+T2{ti)] {a+{h)cT+it2)), 

(45) 

d{a-ih)a-it2))/dh = -iujA{'J~{h)a-{t2)) 

-[ri(ti) + r*(ii)] (a_(ii)a_(i2)), 

(46) 

where ri(ti) and r2(ti) are defined in Eqs. (gD]) and (|¥T|) . 
respectively. One can see that the evolution equations 



of the two-time CF's, Eqs. (|45|) and (l46l) . have the same 
forms as the evolution equations of the single-time expec- 
tation values {a-+{ti)) and (cr_(fi)), Eqs. (gl]) and (^5]) . 
respectively. Hence the QRT holds in this case. 

Case 2: [A, L] = or [B, L{t)] = 0. In this case, using 
Eq. ([25l) . we obtain 

d{a-{ti)<7^{t2))/dti = -iuA{(7-{ti)(T^{t2)) 

-[ri(ti)-Hr;(ii)] (a_(ii)a,(t2)) 

^2T^{h,t2){^z{tl)(T-{t2)), (47) 

d(a,(ii)a_(t2))Mi = - [Ti{ti) + Tlih) + TUti) + r2(ii)] 

x(CT,(ti)CT_(i2)) 

-[ri(<i) + rt(ii)-r*(ti)-r2(ii)] 

x('T-(i2)) 
-2r4(tl,t2)(fT-(<l)f7.(i2)), (48) 

where 

r draih ~ T)e''^^^'--^\ (49) 
Jo 

U{tut2)^ [ ' dT/3{h ~ T)e~'^^^'---l (50) 
Jo 

When we obtain Eq. (|T7| from Eq. p5|) . the last term 
containing r3(ti,t2) in Eq. (H7)l does not vanish since 
[L^ A] ^ and [B,L{t - tz)] ^ 0. Similarly, when we 
obtain Eq. (gH]) from Eq. (gt]), because [L,A] ^ and 
[B, L^(t — ^2)] 7^ 0, the last term containing r4(ti, t2) in 
Eq. (gB]) exists. Thus, compared with single-time equa- 
tions (g5|) and (gi|) . the two-time equations (|T7)) and 
(g5)l have the extra last terms containing r3(ti,t2) and 
r4(ii,i2), respectively. As a result, the QRT does not 
hold in this case. It is also obvious from the sole ap- 
pearance of the individual coefficient of either T^{ti,t2) 
in Eq. (gT]) or r4(ti,i2) in Eq. (gg]) that the finite- 
temperature bath CF's a{ti — t) and /3(ti — t) can not be 
combined into the single effective bath CF acff (ii — t) — 
a{ti — r) + f3{ti — t) of Eq. ([57)1 as in the Hermitian 
coupling operator case. 

Case 3: [A,L^] = or [B,D{t)\ = 0. Using Eq. 
we obtain for this case 

d{a+{ti)a^{t2)) / dti = +iujA{(7+{ti)(T:,{t2)) 

-[r^(ti)+r2(ti)] (a+(ti)a,(i2)) 

-2r4(tl,t2)(a,(ti)f7+(i2)), (51) 

d(a,(ii)a+(t2))Mi - - [ri(ii) + r*(ti) + Tiih) + r2(ti)] 

■x{<J^{ti)(y+{t2)) 

-[ri(ii) + rj(ti)-r*(ii)-r2(ti)] 

x(^T+(t2)) 

-2r3(ii,t2)(a+(ti)a,(i2)). (52) 

Similarly, compared with the single-time evolution equa- 
tions ((42)) and <\AA\ , the two-time evolution equations (jST)) 
and (|5^ have the extra last terms containing T\{ti) and 
r3(ii), respectively. As a result, the QRT also does not 
hold for these two CF's. 
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Case J,: [A, L] ^ 0, [A, L^] ^ and [B, L{t)] ^ , 
[B, L'^{t)] 7^ 0. In this case, by using Eq. (1^ . we obtain 
the following equations 

d{a,{ti)cT,{t2))/dti = - [ri(<i) + T\{ti) + T\{ti) + T2{ti) 

x(CT2(ii)CT^(t2)) 

x{crzit2)) 
+ir3{h,t2){a+{h)a-{t2)) 
+Ari{h,t2){a^{ti)a+{t2)) . (53) 

The evolution equation of the CF {crz{ti)<Tz{t2)), 
Eq. ^S^, is coupled with the evolution equations of the 
CF's (cr+(ii)cr_(t2)) and {(T^{ti)a+{t2)) , which corre- 
spond to the CF's in Case 2 and Case 3, respectively. 
Their evolution equations, obtained from Eq. ([25]), are 



d{a^{ti)a+{t2))/dti = -iujA{cr-{ti)<y+{t2)) 

-[Ti{h)+n(h)] {a.{h)a+{t2)) 
+r3{ti,t2){(7,{h)a,{t2)) , (54) 

d{a+{ti)a-{t2))/dti = iujA{cr+{ti)o-^{t2)) 

-[riih)+T2{h)] {a+{h)a^{t2)) 

+T4{h,t2){<J,{ti)<J,{t2)) . (55) 



Prom Eqs. ([Ml), & and ([55]), it is obvious that the QRT 
also does not hold for CF's {<Tz{ti)<7z{t2)} , (o'+(ii)o'- (^2)) 
and {a-{ti)a+{t2)). 

We may consider any spectral density for which the 
time-convolutionless perturbation theory is still valid to 
characterize the environment, but for simplicity we con- 
sider a spectral density J(w) = J2x I.9aP^('^ — i^a) = 
7?icj(w/A)"~-^ exp(— oj^/A^) with n — 1 (Ohmic), where 
A is the cut-off frequency and 7 is a dimensionless con- 
stant characterizing the interaction strength to the envi- 
ronment. Figure [1] shows the real part of the time evo- 
lution of the system operator CF's (a-+(ti)cr_ (^2)) and 
{crz{ti)az{t2)), as well as the Fourier spectrum S{uj) of 
{a+{ti)a-{t2))- The CF's are obtained in three differ- 
ent cases: the first is in the Markovian case [i.e, taking 
the reservoir CF's a{ti — r) and /3(ti — r) in Eq. ([25]) 
to be 5-correlated in time, or equivalently taking the co- 
efficients of Fi, rj, r2, and t\ to be time- independent 
and equal to their Markovian long-time values and set- 
ting aU T3{ti,t2) and T4{ti,t2) to be zero in Eqs. (l53t . 
([M| and ([55]) ]. the second is in the non-Markovian case 
with a finite cut-off frequency but wrongly directly using 
the QRT method [i.e., the last two terms of Eq. ([^5]) or 
equivalently the terms containing T3{ti,t2) or T4{ti,t2) 
in Eqs. ([53]) . ([54]) and ([55]) being all neglected], and the 
third is in the non-Markovian case with a finite cut- 
off frequency [i.e., using the evolution equation (|25p or 
equivalently Eqs. ([53]) . ([54]) and ([55]) derived in this pa- 
per]. The initial environmental state is in the thermal 
state and the system state in Fig. [Ila)-(c) is arbitrar- 
ily chosen to be |^) = |e) + i \g)^. We can see 




FIG. 1. (Color online) (a) Real part of the time evolution 
and (b) Fourier spectrum S{ij) of the system operator CF 
{a+{ti)a-{t2)), and (c) real part of {az{ti)az{t2)) for three 
different cases: Markovian using the QRT (blue dot-dashed 
line), non-Markovian using the QRT (green dashed line) and 
non-Markovian (red solid line) using the evolution equation 
([25)1 . Other parameters used are uja = 3, (kBT/h) = 1, 
A = 5, 7 = 0.1, and t2 = 1. (d) Fourier spectrum S{ijj) 
of {a+{t-i_)a-{t2)} for a different parameter of 7 = 0.35 and 
for an initial mbced system state. The results of the non- 
Markovian QRT case and the non-Markovian evolution case 
in (c) become indistinguishable when t = ti—t2 is larger than 
1.5. 



that there are considerable differences between the re- 
sults obtained in the three different cases in Fig.[TIa) and 
(b), and more significant differences can be observed in 
Fig- [2c) and (d). The oscillations of the CF's are more 
pronounced in the non-Markovian cases. In Fig. [H^b), 
the coherent peaks of the Fourier spectrum centered at 
cj — ±UJA are higher and the widths are narrower in the 
non-Markovian cases. The CF {(Jz(ti)<Jz{t2)) of the non- 
Markovian evolution equation case (in red solid line) in 
Fig. (He) differs more from the CF's of the other two 
cases (in green dashed line and in blue dot-dashed line) 
in the short-time regime than the CF (cr+(ti)a_(i2)) of 
the non-Markovian evolution equation case (in red solid 
line) in Fig. [Ija) does. This is because as compared to 
the evolution equation of {a+{ti)a-{t2)) of Eq. ([55]) . the 
evolution equation of {(Tz{ti)az{t2)) of Eq. ([55]) has, in 
addition to a term proportional to r4(<i,t2), an extra 
correction term proportional to F3(ti,t2) over its QRT 
counterparts. It is also found that generally the results 
of the non-Markovian QRT case (in green dashed lines) 
approach those of the non-Markovian evolution equation 
case (in red solid lines) more closely than the results of 
the Markovian QRT case (in blue dot-dashed lines) do. 
Similar behaviors are also observed when the tempera- 
ture is increased or when the cut-off frequency A is in- 
creased. The Markovian case can be recovered from the 
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non-Markovian ones in the limit when the cut-off fre- 
quency A — >■ oo, in which the three resuhs coincide. For 
a larger 7 and for an initial mixed system state with the 
values of the off-diagonal density matrix elements being 
a quarter of those of the pure state 1 5') , the peak heights 
of 5(cij) are lower as shown in Fig. [TJd). Furthermore in 
Fig. (rfd), the two coherent peaks are still clearly visible 
in the non-Markovian cases, while the two peaks is barely 
visible in the Markovian case. 

For the present spin-boson model with the system cou- 
pling operator L ^ L\ the self-Hamiltonian of the spin 
does not commute with the system coupling operator, 
i.e, [Hg^L] ^ 0, and the environment coupling oper- 
ator also does not commute with the self-Hamiltonian 
of the environment, i.e., [Hfi^a\\ ^ 0. Thus the exact 
non-Markovian finite-temperature two-time CF's of the 
present spin-boson model are not directly available. But 
in Ref. ^Z] , we evaluated the exact non-Markovian finite- 
temperature two-time CF's of the system operators for 
an exactly solvable pure-dephasing spin-boson model in 
two ways, one by exact direct operator technique without 
any approximation and the other by the derived evolu- 
tion equation (|25|) valid to second order in the system- 
environment interaction Hamiltonian. The perfect agree- 
ment of the results between the non-Markovian evolution 
equation case and the exact operator evaluation case, and 
the significant difference between the non-Markovian evo- 
lution equation case and the case of wrongly applying 
non-Markovian QRT [23. demonstrate clearly the va- 
lidity of the derived evolution equation ((25t . It is thus 



believed that in the weak system-environment coupling 
limit, the finite-temperature CF's calculated using our 
evolution equation that takes into account the nonlocal 
environment memory term, Eq. (|18p , for the present spin- 
boson model would agree more closely with the exact 
results than those in the non-Markoian QRT and Marko- 
vian QRT cases. 



that vanishes in the Markovian case but makes the evo- 
lution equation deviate from the QRT in general cases. 
The derived two-time equations are valid for thermal en- 
vironments at any temperature with Hermitian or non- 
Hermitian coupling operators and for any initially factor- 
ized (separable) system-reservoir state (pure or mixed) as 
long as the assumption of Eq. ([7]) and the approximation 
of the weak system-environment coupling that are used to 
derive the equations apply. In contrast to the evolution 
equations presented in Refs. [T8l-f2^. ^22^ that are applica- 
ble for bosonic environments, Eq. (j22p derived in this pa- 
per can be used to calculate the two-time CF's for a wide 
range of system-environment models with bosonic and/or 
fermionic environments. We have also given conditions 
on which the QRT holds in the weak system-environment 
coupling case and have applied the derived equations to a 
problem of a two-level system (atom) coupled to a finite- 
temperature thermal bosonic environment (electromag- 
netic fields), in which the system coupling operator is 
not Hermitian, L ^ L\ and the evolution equations de- 
rived in Refs. [l"8l - [20j are not applicable. It is easy to 
calculate the two-time CF's using the derived evolution 
equations. Other non-Markovian open quantum system 
models that are not exactly solvable can be proceeded 
in a similar way to obtain the time evolutions of their 
two-time system operator CF's valid to second order in 
the system-environment interaction Hamiltonian. This 
illustrates the practical usage of the evolution equations. 
Therefore, the derived evolution equations that general- 
ize the QRT to the non-Markovian cases will have broad 
applications in many different branches of physics, such 
as quantum optics, statistical physics, chemical physics, 
quantum transport in nanostructure devices and so forth 
when the properties related to the two-time CF's are of 
interests. 
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